/*  sxfbatn2.cpp  by K.Tsuru  */
// function ID 5101  BRADIX
/*************************************************************
SDecimal class
arctan(n/d) by series, n and d are integers.
 arctan(n/d) = (n/d){ 1 - (n^2/d^2)/3 + (n^4/d^4)/5 - .....}
 arctan(n/d) = (n/d) - (n^3/d^3)/3 + (n^5/d^5)/5 - .....
Both series have the same precision.
num < den <= ULONG_MAX/BRADIX = 131071
den > BRADIX is Ok
**************************************************************/
#ifndef SN_H
#include "sn.h"
#endif
static const char* const func = "Batan2";
SDecimal Batan2(ulong num, ulong den){
    RealSize C;
    uint upPrec = 0;
    SDecimal t, res, delta;
    const ulong mt = t.SlOpMaxValue();  //=131071

    if(!num) return res;    // res = 0
    if( (num >= den) || (den > mt) ) t.SetError(t.OUT_OF_RANGE, func, 5101);
    if(num < BRADIX){
        t.SetInt((int)num);     XsDiv(t, den, res);
    } else {
        //When den > num >= BRADIX, the error increases by about two figures.
        //Then it temporally increases the effective figures.
        //Because (SN)*(SN) is not used here, it is not necessary to consider "proper".
        upPrec = 2u;
        C.SetEffFig(C.EffFigures()+upPrec);
        t.SetInt(1);        XsDiv(t, den, res); // res = 1/den
        XsMult(res, num, res);  // res = num/den  < 1.0
    }
    t = res;

//The square of denominator is less than "mt" or not.
    int dsq_lt_mt = ( (double)den*(double)den < (double)mt ) ? 1 : 0;
    ulong k = 1;
    if(dsq_lt_mt){
        ulong dsq = den*den, nsq = num*num;     //square of denominator and numerator
        do {
            XsDiv(t, dsq, t);                   // t /= (den*den);
            if(nsq != 1) XsMult(t, nsq, t);     // t *= num*num
            if( (k += 2) >= mt ){
                t.SetError(t.NOT_CONVERGE, func, -5101);
                break;
            }
            XsDiv(t, k, delta);                 // delta = t / k;
            if(k & 2)   XXSub(res, delta, res); // res -= delta;
            else        XXAdd(res, delta, res); // res += delta;
        } while( delta.Sign() );
    } else {
        int nsq_lt_mt = ( (double)num*(double)num < (double)mt ) ? 1 : 0;
        ulong nsq = nsq_lt_mt ? num*num : 0;
        do {
            XsDiv(t, den, t);   XsDiv(t, den, t);   // t /= (den*den);
            if(num != 1){   // t *= num*num
                if(nsq){
                    XsMult(t, nsq, t);
                } else {
                    XsMult(t, num, t);      XsMult(t, num, t);
                }
            }
            if( (k += 2) >= mt ){
                t.SetError(t.NOT_CONVERGE, func, -5101);
                break;
            }
            XsDiv(t, k, delta);                 // delta = t / k;
            if(k & 2) XXSub(res, delta, res);   // res -= delta;
            else      XXAdd(res, delta, res);   // res += delta;
        } while( delta.Sign() );
    }
    res.upToTerm = (k + 1)/2;
    if(upPrec){
        C.SetEffFig(0);     t = res;    //remove lower figures
        return t;
    }
    return res;
}